Poisson Blending and Mixed-Gradient Blending

This is for the implementing Poisson blending and the blending with the mixed gradient.
I should insert two penguin images on top of the snow in the hiking image.

Load Image

im_background = imresize(im2double(imread('./data/hiking.jpg')), 0.5, 'bilinear');
im_object1 = imresize(im2double(imread('./data/penguin-chick.jpeg')), 0.5, 'bilinear');
im_object2 = imresize(im2double(imread('./data/penguin.jpg')), 0.5, 'bilinear');

Get a mask on the desired part of the source image

obj1_mask = getMask(im_object1);

Draw polygon around source object in clockwise order, q to stop

Specify a location in the target image where it should be blended.

Transform the source image so that indices of pixels in the source and target regions match.

This function was originally given, but I add the downsizing function here for more realistic blending.
In the process of aligning the pixels of source and target, the downsize factor is muliplied to the index of pixels to resize.
In this penguin, I multiply 0.7 times for downsizing.

[im_s1, mask_s1] = alignSource(im_object1, obj1_mask, im_background,0.7);

choose target bottom-center location

The formulation and implementation of poissonBlend is similar to the gradient processing in the toy example.
The difference is as follows.
1. We should consider the gradient values of only "the masked source image"
2. If the certain pixel value is on the boundary of the masked source image, then we should use the target value for computing that of gradients.
3. The pixel values outside the masked region should follow the values from target.

im_blend1 = poissonBlend(im_s1, mask_s1, im_background);

figure(3), hold off, imshow(im_blend1);

The above is the result. The penguin blended much more naturally.

I do the same process with the other penguin.

obj2_mask = getMask(im_object2);

Draw polygon around source object in clockwise order, q to stop

[im_s2, mask_s2] = alignSource(im_object2, obj2_mask, im_blend1,1);

choose target bottom-center location

im_blend2 = poissonBlend(im_s2, mask_s2, im_blend1);

figure(3), hold off, imshow(im_blend2);

imwrite(im_blend1,'Poisson_blend1.png');

imwrite(im_blend2,'Poisson_blend2.png');

The blending is done quite well. But there are some limitations.

Blending is too sensitive to the masking. Because of the imperfect mask, the right side of the first penguin has some white-blocked noise.
The background part inside the mask is blurry. I think it's because the target background is not mixed well in that part.
The second limitation could be resolve through the below mixed gradient.

Mixed-Gradient

mixed gradient blending algorithm has one difference with the poisson blending.
Instead of using the source gradient only as the guide, it use the gradient in source or target with the larget magnitude.
I implemented mixed-gradient algorithm with the same source images.

Penguin1

im_blend_mix1 = mixedBlend(im_s1, mask_s1, im_background);
figure(3), hold off, imshow(im_blend_mix1);

Penguin 2

im_blend_mix2 = mixedBlend(im_s2, mask_s2, im_blend_mix1);

figure(3), hold off, imshow(im_blend_mix2);

The above is the results.
Although the blur inside the source mask appears to be better (goes well with the background of the target image), one serious problem is that the target gradient is projected where it is unnecessary, making the object appear transparent.

imwrite(im_blend_mix1,'Mix_blend1.png');
imwrite(im_blend_mix2,'Mix_blend2.png');

function [output] = poissonBlend(im_s,mask_s,im_background)
    [imh, imw, nn] = size(im_background);
    im2var = zeros(imh, imw);
    im2var(1:imh*imw) = 1:imh*imw;
    
    % mask_1의 x,y좌표값 가져오기 -> A의 ncol 계산하기
    [y,x] = find(mask_s ~= 0);
    [y_s,x_s] = find(mask_s == 0);
    
    % matrix A의 윗부분 만들기 
    nrow = imh*imw;
    ncol = size(x,1)+size(x_s,1);
    A = sparse([], [], [], ncol, nrow, 0);   
    b = sparse(ncol,1);
    output = zeros(imh,imw,nn);
    
    for c = 1:nn
        t = im_background(:,:,c);
        s = im_s(:,:,c);
        e = 0;
 
        % 마스크 안은 source의 gradient를 따르도록 함. 
        % 단, 만약 boundary에 걸쳐있다면, target_j의 값과 source_i의 값 비교 
        for i = 1:size(x,1)
            e = e+1;
            right_term = 0;
            A(e, im2var(y(i),x(i))) = 4;
            
            if mask_s(y(i),x(i)-1) == 0
                A(e, im2var(y(i),x(i)-1)) = 0;
                right_term = right_term + s(y(i),x(i)) - s(y(i),x(i)-1) + t(y(i),x(i)-1);
            else 
                A(e, im2var(y(i),x(i)-1)) = -1;
                right_term = right_term + s(y(i),x(i)) - s(y(i),x(i)-1);
            end
 
            if mask_s(y(i),x(i)+1) == 0
                A(e, im2var(y(i),x(i)+1)) = 0;
                right_term = right_term + s(y(i),x(i)) - s(y(i),x(i)+1) + t(y(i),x(i)+1);
            else 
                A(e, im2var(y(i),x(i)+1)) = -1;
                right_term = right_term + s(y(i),x(i)) - s(y(i),x(i)+1);
            end
            
            if mask_s(y(i)-1,x(i)) == 0
                A(e, im2var(y(i)-1,x(i))) = 0;
                right_term = right_term + s(y(i),x(i)) - s(y(i)-1,x(i)) + t(y(i)-1,x(i));
            else 
                A(e, im2var(y(i)-1,x(i))) = -1;
                right_term = right_term + s(y(i),x(i)) - s(y(i)-1,x(i));
            end
 
            if mask_s(y(i)+1,x(i)) == 0
                A(e, im2var(y(i)+1,x(i))) = 0;
                right_term = right_term + s(y(i),x(i)) - s(y(i)+1,x(i)) + t(y(i)+1,x(i));
            else 
                A(e, im2var(y(i)+1,x(i))) = -1;
                right_term = right_term + s(y(i),x(i)) - s(y(i)+1,x(i));
            end
 
            b(e) = right_term;
        end
    
        %마스크 밖은 target의 pixel을 따르도록 함 
        for i = 1:size(x_s,1)
            e = e+1;
            A(e,im2var(y_s(i),x_s(i))) = 1;
            b(e) = t(y_s(i),x_s(i));
        end
    
        v = A\b;
    
        n=0;
        for j = 1:imw
            for i = 1:imh
                n = n+1;
                output(i,j,c) = v(n);
            end
        end
    
    end
end
 
function [output] = mixedBlend(im_s, mask_s, im_bg)
    [imh, imw, nn] = size(im_bg);
    im2var = zeros(imh, imw);
    im2var(1:imh*imw) = 1:imh*imw;
    
    % mask_1의 x,y좌표값 가져오기 -> A의 ncol 계산하기
    [y,x] = find(mask_s ~= 0);
    [y_s,x_s] = find(mask_s == 0);
    
    % matrix A의 윗부분 만들기 
    nrow = imh*imw;
    ncol = size(x,1) + size(x_s,1);
    A = sparse([], [], [], ncol, nrow, 0);   
    b = sparse(ncol,1);
    output = zeros(imh,imw,nn);
    
    for c = 1:nn
        t = im_bg(:,:,c);
        s = im_s(:,:,c);
        e = 0;
        for i = 1:size(x,1)
            e = e+1;
 
            right_term = 0;
            A(e, im2var(y(i),x(i))) = 4;
 
            if mask_s(y(i),x(i)-1) == 0
                A(e, im2var(y(i),x(i)-1)) = 0;
                right_term = right_term + mix_compare(y(i),x(i),t,s,0,-1) + t(y(i),x(i)-1);
            else 
                A(e, im2var(y(i),x(i)-1)) = -1;
                right_term = right_term + mix_compare(y(i),x(i),t,s,0,-1);
            end
 
            if mask_s(y(i),x(i)+1) == 0
                A(e, im2var(y(i),x(i)+1)) = 0;
                right_term = right_term + mix_compare(y(i),x(i),t,s,0,1) + t(y(i),x(i)+1);
            else 
                A(e, im2var(y(i),x(i)+1)) = -1;
                right_term = right_term + mix_compare(y(i),x(i),t,s,0,1);
            end
            
            if mask_s(y(i)-1,x(i)) == 0
                A(e, im2var(y(i)-1,x(i))) = 0;
                right_term = right_term + mix_compare(y(i),x(i),t,s,-1,0) + t(y(i)-1,x(i));
            else 
                A(e, im2var(y(i)-1,x(i))) = -1;
                right_term = right_term + mix_compare(y(i),x(i),t,s,-1,0);
            end
 
            if mask_s(y(i)+1,x(i)) == 0
                A(e, im2var(y(i)+1,x(i))) = 0;
                right_term = right_term + mix_compare(y(i),x(i),t,s,1,0) + t(y(i)+1,x(i));
            else 
                A(e, im2var(y(i)+1,x(i))) = -1;
                right_term = right_term + mix_compare(y(i),x(i),t,s,1,0);
            end
 
            b(e) = right_term;
        end
    
        for i = 1:size(x_s,1)
            e = e+1;
            A(e,im2var(y_s(i),x_s(i))) = 1;
            b(e) = t(y_s(i),x_s(i));
        end
    
        v = A\b;
    
        n=0;
        for j = 1:imw
            for i = 1:imh
                n = n+1;
                output(i,j,c) = v(n);
            end
        end
    
    end
end
 
function [d] = mix_compare(y,x,t,s,a,b)
 
    if abs(t(y,x)-t(y+a,x+b)) < abs(s(y,x)-s(y+a,x+b))
        d = s(y,x)-s(y+a,x+b);
    else
        d = t(y,x)-t(y+a,x+b);
    end
end
 
function [mask, poly] = getMask(im)
% [mask, poly] = getMask(im)
% Asks user to draw polygon around input image.  Provides binary mask of
% polygon and a chain of all interior boundary points.
 
disp('Draw polygon around source object in clockwise order, q to stop')
figure(1), hold off, imagesc(im), axis image;
sx = [];
sy = [];
while 1
    figure(1)
    [x, y, b] = ginput(1);
    if b=='q'
        break;
    end
    sx(end+1) = x;
    sy(end+1) = y;
    hold on, plot(sx, sy, '*-');
end
 
mask = poly2mask(sx, sy, size(im, 1), size(im, 2));
if nargout>1
    [poly.x, poly.y] = mask2chain(mask);
end
end
 
function [im_s2, mask2] = alignSource(im_s, mask, im_t,a)
% im_s2 = alignSource(im_s, mask, im_t)
% Asks user for bottom-center position and outputs an aligned source image.
 
figure(1), hold off, imagesc(im_s), axis image
figure(2), hold off, imagesc(im_t), axis image
[y, x] = find(mask);
y1 = min(y)-1; y2 = max(y)+1; x1 = min(x)-1; x2 = max(x)+1;
im_s2 = zeros(size(im_t));
disp('choose target bottom-center location')
[tx, ty] = ginput(1);
 
yind = (y1:y2);
yind2 = (yind -max(y))*a + round(ty);
yind2 = round(yind2);
xind = (x1:x2);
xind2 = (xind - round(mean(x)))*a + round(tx);
xind2 = round(xind2);
 
y = (y - max(y))*a + round(ty);
x = (x - round(mean(x)))*a + round(tx);
x = round(x);
y = round(y);
ind = y + (x-1)*size(im_t, 1) ;
mask2 = false(size(im_t, 1), size(im_t, 2));
mask2(ind) = true;
im_s2(yind2, xind2, :) = im_s(yind, xind, :);
im_t(repmat(mask2, [1 1 3])) = im_s2(repmat(mask2, [1 1 3]));
 
figure(1), hold off, imagesc(im_s2), axis image;
figure(2), hold off, imagesc(im_t), axis image;
drawnow;
end